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To develop a sound base on which to design efficient lightguides, 
it is necessary to understand the physics of the propagating modes of 
this new and important transmission medium. To achieve this end, 
we have calculated numerical solutions to the vector form of Max- 
well's equations for lightguides and have obtained all six components 
of the electromagnetic field. We have included materials dispersion 
and calculated Poynting vectors and effective and group indices. 

I. INTRODUCTION 

A possible first step in analyzing production line mcvd (modified 
chemical vapor deposition) lightguide preforms 1 is to calculate the 
effective indices for a propagating wave in a fiber assumed to be 
produced from the preform. A knowledge of the effective indices for 
various wavelengths permits one to calculate the group indices. Using 
these, we can calculate the group delay of the modes. In this way, we 
can ascertain the reproducibility of the preforms and ultimately in- 
crease yields of the products. 

The solution of Maxwell's equations for a radially inhomogeneous 
medium such as is found in lightguides is by no means easy. This is 
because we must deal with fourth-order differential equations with 
variable coefficients. Various approximations can, of course, be 
made. 2 " 8 For example, a popular one is to neglect the gradient of the 
logarithm of the dielectric constant with respect to radial distance 9 
and then employ wkb methods. 10 

These approximations require that discontinuities in index be of 
little importance and that the refractive index vary slowly with dis- 
tance. Solutions of this sort have had a substantial impact on lightguide 
engineering." 
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VERY THICK CLADDING 




Fig. 1— Schematic diagram of a lightguide showing an inhomogeneous core and a very 
thick (essentially infinite) cladding. R or p is the radial coordinate, and <j> is the polar 
angle (p = K<Ji). 

If we cannot or do not wish to limit ourselves to the restrictions 
imposed by these approximations, we have little choice but to keep 
Maxwell's equations intact and resort to numerical techniques. For 
example, Vigants and Schlesinger 12 have analyzed microwave surface 
waves on radially inhomogeneous cylinders using numerical proce- 
dures. In a like manner, Vassel 13 has studied selfoc 14 optical wave- 
guides. 

In this paper, we extend these methods and apply them to the 
analysis of idealized lightguides. In particular, we examine those as- 
pects of the physics of lightguides that are most clearly elucidated by 
the numerical solutions. We shall apply these methods to production- 
line preforms and lightguides in subsequent papers. 

A choice of wavelengths and fiber diameters to study in the present 
paper is necessary. From a purely engineering point of view, a 50- 
micron-diameter fiber and wavelengths of 0.6328, 0.820, 0.850, 1.30 and 
1.55 microns would be essential. This diameter and these wavelengths 



1 1 76 THE BELL SYSTEM TECHNICAL JOURNAL, SEPTEMBER 1 980 



will be applied to all production-line preforms and lightguides analyzed 
in future work. A very large fiber (300 microns) with very many modes 
would be a severe test of our numerical analysis, but would be of little 
practical interest. 

A compromise is clearly necesary. We have chosen to study a 100- 
micron-diameter fiber. This adequately tests our programs and our 
ability to handle many modes. As we are interested in this paper in 
rather uniformly spaced wavelengths, our choice of wavelengths does 
not coincide exactly with those of engineering importance. 

II. THEORY 

We consider a cylindrical lightguide with an inhomogeneous core 
and a thick (essentially infinite) cladding (see Fig. 1). The bound 
propagating modes in the core will be of the form: 



T(R) exp i[wt -M<t>- K a N e Z]. 



(1) 



Here <f> is the polar angle, M the azimuthal mode number, Z the 
longitudinal coordinate, R the radial coordinate, K Q the free-space 
propagation parameter, and N e the effective refractive index of the 
wave. 
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Fig. 2— A possible solution to eq. (2) for Hz. To the left of the line d„ the solution is 
purely oscillatory. To the right of the line, it is the sum of a rising and decaying 
exponential. In this plot, M = 4. 
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Fig. 3— The wavelength variation of refractive index of 13.5-mole-percent germania- 
doped silica. 

Maxwell's equations in cylindrical coordinates can be rigorously 
written as 6,7,15 



^ = p- 1 A(p)r(p). 
dp 



(2) 



Here T(p) is a column vector whose four elements are related to the 
tangential components of the electromagnetic field. These components 
will be continuous through the core-cladding interface. They are as 
follows: 

T 2 = -i p H^Zo (3) 

r 3 = pE* 

I 4 = — llizZio, 

where Z a is the wave impedance of free space and the variable p is 
defined as KM- 

The 4X4 matrix A(p) can be written as: 

(NI/k) - 1 -MN b /k 

pK-M 2 MN e . (4) 

MN e /K p 2 - (M*/k) 

-MN e NI-k 

In this matrix, k is the dielectric constant that varies throughout the 
core of the fiber and remains constant in the cladding. In general, k 
will be expected to change with frequency to. 

It can be shown 1213 that, as p -> », we obtain a form of Bessel's 
equation. For decaying solutions in the cladding, we consequently 
employ modified Bessel functions of the second kind Km(0- There will 
be two independent solution vectors in the cladding. They are 
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Here k(c) is the dielectric constant in the cladding and 

£ - [N 2 e - K(C)] 1/2 PCC 

W M (0 = [N 2 e - K (c)r 4 K M (0 
YAf(0 = £KM0/K M (8, 



(5) 



(6) 



where p cc is the p value at the core-cladding interface. 

To find the solution vectors in the core, we make the following 
change in variable in eq. (2): 

r = p M A. 



We thus obtain the differential equation, 

^ = - [A(p) - MI]A, 
dp p 



(7) 



(8) 



where / is a unit diagonal matrix. 

By requiring that the solution be bounded as p — > o, we obtain the 
two vectors, 13 
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1.476 



Fig. 4— Plot of the function log | G(N e ) \ versus effective index for M = 4 and A = 0.806 
micron. The index profile is parabolic. 
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Fig. 5— A high resolution plot of log | G(N e ) | for the mode marked A in Fig. 4. This 
splitting is due to the presence of the cladding and amounts to a few times 10 . 



'Nl ~ K(0) 

Mk(o) 

MN e 





and 





MN e 

M 

LNl - k(o). 



(9) 



where k(o) is the dielectric constant at the center of the lightguide. 

By using the two solution vectors at the center of the lightguide and 
eq. (8), we can obtain two numerical solutions at the core-cladding 
interface. These must match with the two cladding solutions at the 
same point. This can be written as 

I Tijipcc, N e ) Pj = O, (10) 

where index j specifies a particular vectors and index i a particular 
vector component. Consequently, I\i and I\| are core solutions and 
r«,3 and r,, 4 are cladding solutions. Pi and P 2 are the weights applied 
to the core solutions, and P 3 and P 4 are the weights applied to the 
cladding solutions. The characteristic equation from which we can 
determine N e can then be written as 



G(N e ) = Det[r,> cc , N e )] = o. 



(11) 



This is an implicit, transcendental equation, and to find its solution we 
plot G(N e ) as a function of N e . The zeros will give us the allowed values 
of effective index. 
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III. COMPUTING PROCEDURE 

We obtain numerical solutions to eq. (8) by means of an optimized 
second-order Runge-Kutta procedure. 16 Double-precision arithmetic is 
employed and needed. To minimize computing time, we wrote our own 
programs especially tailored to the above problem. 

There will be regions in the lightguide where the solution to eq. (2) 
will oscillate and other regions where the solution will grow or decay 
exponentially. As we are limited by our computer to numbers between 
10 -38 to 10 +38 , this will limit the wavelength and size of the fiber we can 
study. However, as will be seen, realistic fibers and wavelengths can 
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Fig. 6 — A plot of the absolute value of the splitting of the modal pairs vs mode 
number Q. The dramatic change in splitting for Q = 22 and 23 is due to the cladding. 
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Fig. 7— Plot of the function log | G{N e ) | vs effective index for M = 4 and X = 1.540 
microns. Note the increase in spacing between the effective indices over the data on 
Fig. 4. 



m 

be handled. As an illustration, Fig. 2 shows a possible solution to eq. 
(2) for H z . To the left of the vertical line d„ the behavior is clearly 
oscillatory. To the right of the line, we have the sum of a rising and 
decaying exponential. The decaying term dominates nearest the line, 
while the rising term dominates further away. The line d s locates the 
turning point, and it should not be confused with the core-cladding 
interface. It is simply that radial position where the character of the 
solution changes from oscillatory to exponential. For loosely bound 
modes, it will be very near the core-cladding interface. For tightly 
bound modes, it will be near the center of the lightguide. 

As mentioned before, the solution in the cladding can be written in 
terms of modified Bessel functions of the second kind. We obtain these 
functions by a double-precision routine written by Gatto and Seery. 17 
Again, because of the range of numbers the computer can handle, 
there is some limitation on the size of the argument the program can 
evaluate. However, realistic wavelengths and fiber diameters are not 
precluded. 

The function G(N e ) that we must plot to obtain the N e 's (i.e., 
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effective indices) has a large dynamic range. We are primarily inter- 
ested in the zeros of this function. The large values can be compressed 
by plotting the logarithms of the absolute values of G(N e ). It should be 
borne in mind that exact zeros for | G(N e ) | are not very likely, as G(N e ) 
is the result of a numerical calculation. What we find are approaches 
to zero, and this causes no difficulty as these approaches are very 
evident. 

IV. EFFECTIVE INDICES FOR AN IDEALIZED FIBER 

We study an idealized lightguide with a 13.5-mole percent germania- 
doped silica core and a pure silica cladding. The core is assumed to be 
perfectly cylindrical with a diameter of 100 microns. The cladding will 
be considered essentially infinite. The doping profile is such that the 
refractive index is a parabolic function of distance R. Thus, 

N 2 -Ni 



N=Ni + 



Re 



& 



(12) 



where N\ is the refractive index at the center of the core, N2 the 
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Fig. 8 — A family of dispersion curves for the various modes. The curves for the core 
and cladding are also drawn. Modes that touch the cladding curve are no longer bound 
to the core. 
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Fig. 9— A plot of effective index vs mode number Q. This highly linear behavior is 
only for purely parabolic index profile. 

refractive index in the cladding, and R cc the radius of the core. Both 
Ni and N 2 will be functions of wavelength. 

Fleming 1819 has accurately determined the wave-length variation of 
the refractive index of silica and some doped silicas. The data for 13.5- 
mole percent germania-doped silica is shown in Fig. 3. As expected, 
the refractive index falls off slowly with wavelength. 

It is necessary to have a formula for the wave-length dependence of 
index data such as is shown in Fig. 3. One which is linear in its 
constants is 20 

A/, = Co + CiX 2 + C 2 A 4 + C 3 /(A 2 - 0.035) 

+ C 4 /(A 2 - 0.035) 2 + C 5 /(A 2 - 0.035) 3 , (13) 

where A is in microns and i = 1, 2. We find that this formula fits the 
data of Fig. 3 with a correlation coefficient 21 of 0.9999. Table I lists the 
coefficients for eq. (1) for both 13-mole-percent germania-doped silica 
and silica. 

There is a physical interpretation of eq. (13). The first terms are a 
Taylor series corresponding to the influence of a far infrared band. 
The last terms are a Laurent series corresponding to an absorption 
band near A = 1870A. We find that this formula is also useful in fitting 
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Table I — Coefficients for eq. 
(1 3) for core and cladding 

Core' Cladding* 



C„ = 1.4706868 
Ci = -0.0026870 
C 2 = -0.0000356 
C 3 = 0.0035756 
C 4 - -0.0000828 
Cs = 0.0000018 



1.4508554 
-0.0031268 
-0.0000381 

0.0030270 
-0.0000779 

0.0000018 



* The deviation of any index from 
the fitted curve is less than 10~ 7 . 



the effective index data for the various modes. For this application, 
the correlation coefficient is equally good. 

Figure 4 shows a plot of the log | G(N e ) | vs effective index for our 
idealized fiber with M = 4 and A = 0.806 micron. We notice that there 
are 23 clear approaches to zero, and these are the allowed values of N e . 
The spacing between the N e values is essentially uniform, and they 
span a range from the core center index (1.47453235) to the cladding 
index (1.4535429). 




1.464 



Fig 10— Plot of the function log | G(N e ) | vs effective index for M = 12 and X = 1.540. 
A comparison with Fig. 7 shows a gap in the allowed effective index values. As M gets 
even larger, the gap will increase further. 
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Table II — Listing of 

the wavelengths* 

used in Figs. 8 and 9 



J 


Xt (microns) 


1 


0.6040 


2 


0.6520 


3 


0.7000 


4 


0.7600 


5 


0.8060 


6 


0.8480 


7 


0.9020 


8 


0.9500 


9 


1.0040 


10 


1.0400 


11 


1.1000 


12 


1.1560 


13 


1.2020 


14 


1.2520 


15 


1.2980 


16 


1.3520 


17 


1.4000 


18 


1.4556 


19 


1.5040 


20 


1.5400 


21 


1.6000 



* Rounded from double 
precision values. 

If we examine the approaches to zero under higher resolution than 
is shown in Fig. 4, we find that, except for the mode marked D, the 
allowed values of N e are doublets. Figure 5 shows a high resolution 
plot of log|G(iV e ) | for the mode marked A. The two values of the 
effective indices are 

Nl = 1.45384520 

N 2 e = 1.45384788. 

The difference between these two values is a few times 10 -6 . This 
splitting is due to the presence of the cladding. For example, a mode 
whose electromagnetic field has decayed before reaching the cladding 
(such as B in Fig. 4) has a splitting of a few times 10~ 7 . Figure 6 is a 
plot of the absolute value of this splitting vs mode number Q. Q counts 
the number of pairs of modes where Q = 1 corresponds to the 
nondegenerate mode D. 

From data of the sort shown in Fig. 4 except that they are taken at 
different wavelengths, it is possible to draw a family of dispersion 
curves for the effective indices. Figure 7 shows a curve analogous to 
Fig. 4 except for A = 1.540 microns. The only obvious change is that 
the spacing between the effective indices has increased. We have 
collected data at 21 wavelengths, and these are listed in Table II. The 
family of dispersion curves that results from this data is plotted in Fig. 
8. Included in the plot are the dispersion curves for the core and 
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cladding material. Note that the resolution of Fig. 8 is not high enough 
to show the splitting of the degenerate modes. 

There is some similarity in shape between the effective index curves 
and the core and cladding curves. We also clearly see the increase in 
spacing between the effective index curves for longer wavelengths. 
Some modes are seen to approach the curve for the cladding. When 
this happens, such modes are no longer bound to the core. Plots like 
Fig. 8 display rather well the characteristics of the lightguide. Another 
sort of plot of the data of Fig. 8 is also informative. Rather than 
plotting N e vs wavelength for a given mode, we can plot N e vs mode 
number Q for a given wavelength. Such a plot is shown in Fig. 9, and 
the wavelengths are those of Table II. The most striking feature is 
that N e for a given wavelength is nearly a linear function of Q. 

Another important aspect of the basic physics of lightguides is 
shown in Fig. 10. This plot is identical with Fig. 7, except that M now 
equals 12. We note that there is a gap in the allowed effective index 
values. Such a gap also exists in the data of Fig. 7, but it is not nearly 
so evident. As M increases further, the gap will get larger, until finally 
all modes are prohibited. 




160 240 

RADIAL DISTANCE 



Fig. 11— The H R wavefunction for M = 4 and A = 0.806 micron. The effective index 
is 1.4538452. Only the field in the core is plotted. The radial distance is in terms of p. 
Note that, at the core cladding interface, the field has a substantial value. 
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Fig 12 — A plot of all six components of the electromagnetic field. Again M = 4, 
a = 0.806, and N e = 1.4538452. 



V. WAVEFUNCTIONS FOR THE ELECTROMAGNETIC FIELD 

Knowing the effective index for a mode, we can generate the wave- 
functions for the electromagnetic field by means of eqs. (3) and (10). 
We again consider the pair of modes labeled A in Fig. 4. Figure 11 
shows the H R wave-function for the lower of the two effective indices. 
The radial distance is in terms of the variable p. The function plotted 
corresponds to the field in the core. In the cladding, which is to the 
right of the dotted line at p = 380, the field decays in a roughly 
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exponential manner. Since Hr is not tangential to the core-cladding 
interface, it will, of course, show a discontinuity at this boundary. 

We note a number of features in Fig. 11. First, the field oscillates 
with decreasing frequency as a function of p. Second, the field has a 
substantial, nonzero value at the core-cladding interface. This means 
that, for this mode, appreciable energy is carried in the cladding. 
Finally, we note that the oscillations for this component do not show 
a continual decrease in amplitude with p. There is an initial decrease 
to a constant level and, just at the core-cladding interface, a slight 
increase. All six components are shown in Fig. 12. Only E z and H z 
decay continuously with p. 

From a knowledge of the field functions, the Poynting vector 22 may 
be calculated. Thus: 



P = -{E R Ht-E„H%). 



(14) 



This vector, which gives us the power density in the guide, is plotted 
vs p in Fig. 13 for the field functions of Fig. 12. We note a decrease in 
power from the guide center to a relatively constant level and then an 
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Fig. 13 — Poynting vector for the field functions shown in Fig. 12. The vertical scale is 
relative. 
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Fig. 14— Poynting vectors for the pair of modes designated C in Fig. 4. We note that 
one mode has three peaks and the other two. Q = 3 for this pair of modes. 



increase at the core-cladding interface. We also note an increase in the 
width of the power peaks as the cladding is approached. 

We may use the Poynting vector to classify the more tightly bound 
modes. By tightly bound, we mean those modes that are little affected 
by the cladding. Modes with Q = 1 through Q = 20 (see Fig. 6) are 
examples of these. To illustrate this classification scheme, consider the 
pair C, (Q = 3) in Fig. 4. The effective indices for these are: 

Nl = 1.47106909 

N 2 e = 1.47106947. 

The Poynting vectors are plotted in Fig. 14 for this pair of modes. The 
core-cladding interface is well to the right of these functions and is not 
shown. Obviously, little energy is being carried in the cladding. We see 
that the mode with the highest efffective index has two peaks and the 
mode with the lowest has three peaks. We designate the latter mode 
as HE 43 and the former EH42. A consistent classification 23 for the two 
modes of each pair is HE m ,q and EH m ,q-\. No EH mode is associated 
with an HE M ,\ mode. 
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VI. CALCULATION OF GROUP INDICES 

A parameter of some importance in lightguide engineering is the 
group index. This will be expected to vary with wavelength and mode. 
Because many modes will normally be excited and because a source 
has some spectral width, there will be a spread in group indices hN g . 
For a fiber of length L, this results in a spread in propagation time At 
given by 



At = ^AN e . 

The group index is related to the effective index by 

N e = N e -\(dN e /d\). 



(15) 



(16) 



This formula may be applied to bulk glass as well, in which case N e is 
simply the usual index of refraction. To apply eq. (16), we must 
calculate dN e /d\ for numerical data, and this could be troublesome. 
We represent the data on both bulk glass and the effective index of 
Fig. 8 by eq. (13). The coefficients are determined by a least-squares 
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Fig. 15 — Group index as a function of wavelength for the bulk germania-doped silica 
glass whose index of refraction dependence is given in Fig. 3. Note that, near 1.4 microns, 
there is very little change of group index with wavelength. 
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Fig. 16 — The family of group index curves for the effective index data of Fig. 8. Note 
the similarity with the curve of Fig. 15. The curves deviating from the bundle are 
tending to become unbound. 

analysis, and this causes a degree of smoothing. The derivative dN e / 
d\ is then easily evaluated. 

Figure 15 shows the group index as a function of wavelength for the 
bulk germania-doped glass data shown in Fig. 3. Of course, since we 
are dealing with bulk glass only, there will be no modal dependence. In 
the neighborhood of 1.4 microns, there is little change in group index 
with wavelength. Thus, at this wavelength a rather broad source would 
transmit data through the bulk material efficiently. 

In Fig. 16, we show the family of group index curves for the data of 
Fig. 7. Each curve corresponds to a particular mode. We note that the 
modes form a bundle with a shape nearly identical to that of the bulk 
glass data of Fig. 15. Because of the steepness of the function near 0.8 
micron, one tends to get a false impression of the relative width of the 
bundle. In fact, one finds that it is only slightly wider at the longer 
wavelengths. The bundle has curves deviating from it, and these are 
modes that are tending to become unbound. From Fig. 16 we estimate 
a spread in group index of about 6 X 10" 4 at 1.4 microns. The time 
spread calculated by eq. (15) is thus 2 X 10" 9 s/km for a broad source. 
The bundle shows narrowing at about 1 micron, and this might be a 
desirable wavelength to operate a narrower source." 
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VII. GRAPHICAL DISPLAY OF GROUP INDEX 

It is useful at times to have a graphical display of group index as a 
function of effective index. This may be helpful in understanding the 
effect of distortions in the index profile on differential modal delay. 
Differentiating eq. (9) yields 

(17. 



— - d\ + 
SK dN e 



or 



dN e 
d\ 



dG 
~d~\, 



dN e = 

I dG 

' dNe 



Thus we obtain the equation 



N e = N e + X[ — 



dG /dG 



d\ dN e 



(18) 



(19) 



This function is plotted in Fig. 17 for A = 1.39 microns, M = 4, and, as 
before, a parabolic index profile. From Fig. 17, we see that N e is made 
up of a family of "dispersion-like curves," and there is one curve for 
each pair of nearly degenerate modes. For the modes that are little 
influenced by the cladding, the curves are very similar. The only 
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Fig. 17 — Plot of eq. (17) relating group index to effective index for A = 1.39 microns. 
The dots correspond to the allowed values of effective index and were taken from plots 
of the sort shown in Figs. 4 and 7. 
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Fig. 18 — Plot of eq. (17) relating group index to effective index for A = 0.803 micron. 
The dots correspond to the allowed values of effective index and were taken from plots 
of the sort shown in Figs. 4 and 7. 



difference is a downward shift as N e gets smaller. For those modes that 
are influenced by the cladding, i.e., those on the extreme left of the 
plot, the curves tend to reverse. 

The dots on the curves are the values of effective index for each 
mode. These data were obtained from plots of the sort described 
earlier (see Figs. 4 and 7). For each mode, we can get a group index. 
Figure 18 shows a plot for the same fiber except at a shorter wave- 
length. In this case, we have A = 0.803 micron, and there are many 
more modes. The curves are much steeper and, as before, they change 
in character as the cladding is approached. 

To accurately calculate the group indices at a particular wavelength, 
one can use either eq. (16) or eq. (19). If we use eq. (16), we must very 
precisely determine the effective indices for the modes of interest at 
three equally spaced wavelengths, for example, 0.81, 0.82, and 0.83 
microns. The required derivative is then easily evaluated. This is a 
simple procedure, but has the disadvantage that precise eigenvalues 
are needed at a number of wavelengths. This is costly in computer 
time. On the other hand, if we use eq. (19), then we only need to 
accurately determine the effective indices at one wavelength. The 
actual calculation of the group index just requires a trivial incrase in 
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computer time over that needed to accurately determine the effective 
indices at a single wavelength. 23 The disadvantage of this procedure is 
that the ratio of two derivatives need be calculated, and this could be 
less accurate. 

VIM. CONCLUSIONS 

Using numerical methods, we have obtained solutions to the vector 
form of Maxwell's equations for idealized lightguides. The results 
obtained clarify the physics and suggest that application of these 
methods to production-line preforms and lightguides will be useful. In 
particular, displays of the effective index data such as shown in Fig. 4 
and the displays of group index data such as shown in Fig. 18 may be 
valuable. In addition to the work described in this paper, we have 
tested the mathematical methods for step-index, single-mode fibers. 
Agreement with known results is excellent. 
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